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Burial of the polar magnetic field of an accreting neutron 
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ABSTRACT 

The hydromagnetic structure of a neutron star accreting symmetrically at both mag- 
netic poles is calculated as a function of accreted mass, Ma, and polar cap radius, 
starting from a centered magnetic dipole and evolving through a quasistatic sequence 
of two-dimensional, Grad-Shafranov equilibria. The calculation is the first to track 
fully the growth of high-order magnetic multipoles, due to equatorward hydromag- 
netic spreading, while simultaneously preserving flux freezing and a self-consistent 
mass-flux distribution. Equilibria are constructed numerically by an iterative scheme 
and analytically by Green functions. Two key results are obtained, with implications 
for recycled pulsars, (i) The mass required to signiflcantly reduce the magnetic dipole 
moment, 10~^ Mq, greatly exceeds previous estimates (~ 10~^'^Mq), which ignored 
the conflning stress exerted by the compressed equatorial magnetic fleld. (ii) Magnetic 
bubbles, disconnected from the stellar surface, form in the later stages of accretion 
(Ma > IO-^Mq). 
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1 INTRODUCTION 

Observations of low-field neutron stars in binary systems 
containing white-dwarf and supergiant companions, with 
a history of disc-fed and wind-fed accretion respectively, 
suggest that the magnetic dipole moment |m| of a neu- 
tron star decreases mo notonical ly with accreted mass, Ma 
jT^Mn_&vai^«i^Ieuyol 1986; van den Heuvel &: Bitzarakil 
I^Qsi rsee IWiierj|l997r for a dissenting view). Several mech- 
anisms have been proposed to explain why |m| is re- 
duced: (i) accelerated Ohmic decay, where the electrical 
conducti vity of the crust is lowered by accretion-induced 
heating iUrpin fc GeDPerti ri995: Uroin & Kononkov 1997); 
(ii) interactions between superfluid neutron vortices and 
superconducting magnetic fluxoids in the stella r interior 
jMusUmov fc Tsvganlil985HSrinivasan et al.lll990t) : and (ifi) 
magnetic screening or burial, where the currents generat- 
ing the natal magnetic field are partial ly neutralized by 
accretion-induced curre nts in the crust iBlondin fc Freesa 
ll986l:lArons fc Leall980l). For a critical review of these mech- 
anisms. see iMelc^o^^PhinnevI i200ll) . 

In this paper, we study the mechanism of magnetic 
burial in detail. In the early stages of accretion (Ma < 
1O~^°M0), accreted matter accumulates in a column at 
the polar cap, minimally distorting the magnetic field. The 
mass-fiux distribution in this regime has been calculated by 
Grad-Shafranov methods, with the prediction that |m| is re- 



duced by ^ 1 per cent for Ma ~ 10 ^°Mr7» j|Hanieurv et alJ 
ll983HBrown fc Bildstenlll998l:lLitwin et alJbOOlT) . We show 
that these calculations overestimate the amount of screen- 
ing; in fact. Ma > 10~^ is required to reduce |m| by 10 per 
cent when the confining stress of the compressed, equatorial 
magnetic field is modelled faithfully. In this regime, inac- 
cessible to previous analyses due to numerical breakdown 
jHameurv et al.lll983l:lBrown fc Bildstenlll998l:lLitwin et all 
l20^F ~the latitudinal pressure gradient at the base of the 
polar column forces the polar magnetic field to buckle, and 
the accreted material spreads equatorward toge ther with 
frozen- in magnetic fiux iMelatos fc Phinnevl200ir) . We com- 
pute the structure of the highly distorted magnetic field, and 
hence |m|, as a function of Ma.. 

A key advance in the present work is that the mass- 
flux distribution in each equilibrium state is self-consistent; 
our equilibria are generated by a continuous deformation 
of the fiux surfaces of the initial magnetic field (say, a 
dipole), in a manner which preserves fiux-freezing. This 
is not true of previous calculations, where the mass-fiux 
distribution is unconstrained rel ative to the initia l state 
( Brow n fc Bildsten 1998; Hamcu rv et al 1IT983: Litwin et alJ 
120011: iMeTatos fc Phinnevll200lh . However, several other im- 
portant effects are not included to keep the problem man- 
ageable, (i) Ohmic dissipation is neglected, even though 
the diffusion and accretion time-scales are comparable for 
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the smallest magnetic structures predicted by the theory 
JCumming et al.ll200ll: ISrown fc Bildstenlll998() . (ii) We do 
not inyestigate the stability of the hydromagnetic equi- 
libria we compute; sharp magnetic-field gradients are po- 
tentially d isrupted by Ray leigh- Taylor and interchange in- 
stabilities llCumming et al.ll2n0ll iMelatos fc Phinnevll2n0ll: 



lBhattacharya^l999^ . (iii) We treat the neutron star as a hard 
surface; subsidence of accret ed material, and inc orporation 
into the crust, are neg lected (jBhattacharvaill999h . 

The paper is structured as follows. In Section 2, we 
introduce the theoretical framework for calculating the self- 
consistent hydromagnetic equilibrium state of an accreting 
neutron star. Analytic and numerical methods of solution 
are giyen in Section 3. The properties of the equilibria are 
inyestigated in Section 4, |m| is computed as a function of 
Ma and the radius of the polar cap, and the noyel feature 
of magnetic bubbles is explored. The limitations of our re- 
sults, with respect to time-dependent processes like hydro- 
magnetic instabilities and ohmic dissipation, are assessed in 
Section 5. 



2 THEORY OF EQUILIBRIA 

2.1 Hydromagnetic force balance 

The equati ons of non-ideal magne tohydrodynamics (MHD) 
in SI units iBernstein et alJll958r) comprise the equation of 
mass conservation. 



dp 
dt 

the equation of motion, 
D'v 

P-^ + p(v- V)v = 



+ V ■ (pv) = , 



(1) 



Vp+— (VxB)xB, (2) 

MO 



and the induction equation (minus the displacement cur- 
rent). 



^-Vx(vxB) = — V'B, 



dt 



^J,oa 



(3) 



supplemented by V • B = and an adiabatic or isothermal 
equation of state, d{pp^^)/dt — 0. In this notation, B, p, 
p, 4>, V and (T represent the magnetic field, mass density, 
kinetic pressure, gravitational potential, plasma bulk veloc- 
ity and electri cal conductivity respectively. Elastic stre sses 
are neglected (lRomanill l990HMelatos fc Phinnevir200ll) . as 
is the Hall effect i Geppert fc Rheinhard1^l2002^ . 

In the magnetostatic limit, defined by v = and d/dt = 
0, the equation of motion reduces to 

Vp + pV(?!)- — (V X B) X B = 0. (4) 

MO 

The local Alfven time-scale, ta ~ L/va ^ 4 x 10~^s for 
(L < 50 m.), is much shorter than the accretion time, 
Ta ~ lO'^ yr. Equations and Q are also satisfied identi- 
cally (in the ideal-MHD limit o — > oo) and drop out of the 
problem. To preserve the information encoded in Q and 
iPJ, we must impose an auxiliary constraint on the mass- 
flux distribution of the final state in order to connect it with 
the initial state and uniquely specify the problem. The con- 
straint expresses the fact that material cannot flow across 
magnetic flux surfaces in the limit cr ^ oo. We delay con- 
sideration of ohmic dissipation, where magnetic flux diffuses 



through the fiuid at short length-scales via (^J, to a future 
paper. 

We define spherical polar coordinates (r, 6, (j)) such that 
^ = defines the symmetry axis of the pre-accretion mag- 
netic field. For an axisymmetric configuration, there exists 
a scalar flux function ip{r, 6) that generates B via 



B 



X e^. 



(5) 



The toroidal component is zero at all times, if the ac- 
cretion process is axisymmetric and = in the initial 
accretion state. Upon substituting jnj into I^J, we obtain 



with 



poi'^ sin^ ( 



sin 6 d 



(6) 



(7) 



We can then resolve Q into components parallel and per- 
pendicular to the magnetic field: 



pVifi + Vp = , 



pV4> + Vp + (A^i/') v^A = . 



(8) 



(9) 



In this paper, we assume the accreted material forms 
an isothermal atmosphere, with p = c^p, where Ca denotes 
the isothermal sound speed. [The force equation for a gen- 
eral equation of state p = p(p) is given in Appendix A of 
^ouschovias (1974).] The gravitational potential (j), deter- 
mined by Poisson's equation, \7^cj> — -iTiGp, is the sum of 
contributions from the accreted material (Ma) and the un- 
derlying neutron star (M,), with Ma <^ M*. As the hydro- 
magnetic length-scale |B|/|VB[ is much smaller than the hy- 
drostatic length-scale |p|/|pV(?!)| (verified a posteriori), V0 is 
approximately constant near the stellar surface for our pur- 
poses, i.e. 



GM,r/Rl 



(10) 



where M, and _R, are the mass and radius of the neutron 
star. 

We use the method of characteristics to solve ||SJ assum- 
ing the gravitational field is radial (Ma ^ M»). The r com- 
ponent reads pr -\- {/^^ip) / c^ipr ~ —p/c^4>r and the 9 com- 
ponent reads pe = {/^"^^l)) / c^tpg, where subscripts indicate 
differentiation. Together these become: pr + {'ipr / '4'e) Pe = 
—p/Cs4>r- The characteristic equation is: dr = {ipg/ipr)d9 = 
—dp/{(f)p). This is solved to yield the two characteristic 
curves: logp + ^/Cg — Ci and ip — C2- Thus the charac- 
teristic solution is logp-|- 4>/c^ — f{ip) or, equivalently. 



p = F(V')exp[-(0-0o)M'] 



(11) 



where F{ip) = exp[f{tp)] is an arbitrary positive function 
to be specified and ~ GM,/R, is a reference potential. 
This is just the usual barometric formula with a different 
base pressure F{ip) for each field line. Note that VF is par- 
allel to VV", so tp and hence F are constant along a field 
line, and one has VF = F'{ip)'V'ip. Substituting into (|SJ, we 
obtain a second order, non-linear, elliptic partial differential 
equation, the Grad-Shafranov equation, for tp: 



-F'(V)exph(,^-0o)/ce'l. 



(12) 
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Equation 1121 can be understood as follows. The quan- 
tity F is a, function of r and 8 through tp at hydrostatic 
equilibrium, expressing the fact that magnetic forces act 
only perpendicular to field lines, while pressure gradients 
balance gravity along field lines. If ^{r,9) is given, and if 
matter is distributed between field lines so that the forces 
parallel to field lines are in exact balance, then forces per- 
pendicular to the field lines are brought into balance by the 
appropriate current density fJ-g^V x B. 

2.2 Magnetic flux freezing and mass-flux ratio 

Many authors guess F{il}) when modelling various sys- 
tems, e.g. st ructures in the solar coro na, like prominences 
and arcades l|Duriggv ll l953l : | L9w| |l980l), and accreting com- 
pact objects (Hamcurv ct al. 1983; Uchida ^ ^owl ligsH: 
iBrown fc B'iidstea.l99a:.,Melatos fc Phinnev-200lh . As F{tp) 
does not change in passing from the initial (pre-accretion) to 
the final (post-accretion) state in the ideal-MHD limit, the 
guessed F{tl}) conflicts with the initial F{il)) except under 
very special circumstances. In this paper, we adopt a self- 
consistent approach whereby we calculate -/^("i/") explicitly by 
demanding that the mass-fiux distribution of the final state 
equals that of the initial state, plus the accreted material, 
as described below. 

Define a local coordinate system (s, t), with unit vectors 
Bs — B/|B| and et — Vtp/\Vtp\ parallel and perpendicular 
to the magnetic field respectively. The amount of matter 
between two infinitesimally separated flux surfaces i/> and 
-i/j -|- dtp is dM = 2tt dsdtpr sin 0, where dsdt = \esds x 
etdt\ = dsd'i/'/l Vi/'l is an inflnitesimal area element, and C 
is the curve tp[r{s),9{s)] — tp. Hence we can write 

dM 



■2n I dsp[r{s),e(s)]rsine\Vtlj\' 
Ic 



(13) 



where dM/dip is the mass-flux ratio. Upon substituting lllH 
into 1131 1. we arrive at 



2tv dip 



ds r sin 0| Vt/)] 



(14) 



which is to be solved simultaneously with 1121 for tp{r,0) 
given dM/dip. For problems involving accretion, the mass- 
flux constraint is not conservative; dM/dtp in the final state 
equals dM/dtp in the initial state plus the mass-fiux distri- 
bution of the accreted material. 

In disc-fed accretion, mass accretes onto polar magnetic 
field lines that close beyond the inner edge of the accretion 
disc, located at a radius 



Ra. 



270 



M 



-2/7 



( |m| 



V 1020 Tm'^ 



4/7 



(15) 



10-^MQyr-\ 

in the equatorial plane ijOhosh fc Lamb| llfl 
iBasko fc SunvaevI Il976ll . The fiux surface that closes 
at Ra. is related to the flux surface at the stellar equator, 
tp* = tp{Rt.,TT/2), by tps. = tptRt/Ra, for a dipole. We 
do not model the mechanism by which plasma enters the 
magnetosphere (e.g. Rayleigh- Taylor and Kelvin-Helmholtz 
instabilities at 7? « -Ra), which sets the form of dM/dtp in 
reality, as this is an unsolved problem iBasko fc SunvaevI 



Il976l: lArons et all Il984l) . Instead, we assume that the 
accreted mass is distributed nearly uniformly within the 
polar fiux tube < tp < tpa_, and that leakage onto flux 
surfaces tpa, < tp < tpt, is minimal. A step change in 
dM/dtp at tpa leads to numerical problems, so we ap- 
proximate the mass distribution over one hemisphere by 
M{tp) = Ma(l - e-'''/''''')/2(l - e-'f"^'''-) . We have checked 
that the solution of 1121 and 1141 is not sensitive to the 
exact functional form of dM/dtp. Finally, we assume that 
the accreted material does not transport any magne tic fiux, 
e.g. from the accretion disc (cf. lUchida fc Lowlll98lll . If the 
magnetic dipole moment is less t han « 10^^ T m'^, we have 
tps. ~ tpt, dCheng fc Zhanj|l99d) and the above functional 
form of M{tp) is inadequate. This occurs at the latest stages 
of accretion (Ma > O.IMq), outside the regime modelled in 
this paper. 



2.3 Initial and boundary conditions 

In this paper, we investigate the distortion of an initially 
dipolar magnetic field. 



tpi{r,9) = tptR-tr ^ sin^ 



(16) 



with tpf = Bt,R^/2, where _B« is the polar magnetic field 
before accretion. Given dM/dtp as a function of M^, we solve 
1121 and 1141 subject to the Dirichlet boundary conditions 



tp{R,,e) =tp,sine and lim tp{r,e) = 0. 



(17) 



the field is dipolar. Far from the star, one has 
as for any localised, static current distribution. 



At the surface, i.e. the crystalline layers of density < 4 x 10^* 
kg m~'^ 
tp (X r 

The approximation that the surface field remains dipo- 
lar at all times is valid provided that Ma is small com- 
pared to M^,, for then the footpoints of the magnetic field 
lines are anchored to the highly conducting, high-inertia 
interior of the star. The surface field may be generated 
deep in t he neutron star core or by a dynamo in the in- 
ner crust I^ Thompson fc Duncanll9 93VKon enkov fc Gepper3 
\200X)- This line- tying boundary condition is a feature of 
models of magnetic loops in the solar corona jLowlll98d: 

and earlier work on neutron 



star accretion juchida & Lowl 119811 


Hameurv et al.l Il983l: 


iBrown fc Bildsten 1998; Litwin et alj 


200l|l. A drawback of 



preventing the accreted matter from sinking is that unreal- 
istically high densities (> 4 x lO^** kg m""^) are produced 
locally, at the base of the column, for Ma > 10~*Mq. Re- 
cent modelling of the magnetic field beneath the surface of 
an accreting neutron star, based on J^J with a velocity dis- 
tribution for the superfluid assumed, illustrates the effects 
of submergence and subsequent incorporation of accreted 
matter into the crust JChoudhuri fc Konaj2o'o3) . ICummi'n3 
(2002) discusses a vacuum (like in this paper) and a screened 
boundary condition at the surface. 

In the numerical calculations presented in Section |^ 
two extra grid-related boundaries are introduced: the outer 
radius of the grid, at r = _Rni, and the lines 9 = and 
9 = ±7r/2, arising when the grid is restricted to one quad- 
rant or hemisphere. We choose Rui large enough to in- 
clude the layer where the accretion-induced screening cur- 
rents lie, i.e., above the greater of the hydrostatic {(^/g) 
and Alfven (|B|/|VB|) scale heights. In practice, this is 
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achieved by increasing i?m until the dipole moment of the 
solution varies by less than 0.1%. At r = Rm, the mag- 
netic field is taken to be radial, with dtp/dr{Rni,6) = 0, 
i.e. a free boundary. (It is outside the scope of this paper 
to model the disc-maeneto sphere interface in detail; (see 
iRastatter Sz Schindleil[l999l) . Another possible way to treat 
the boundary condition at r = Rm is to set 'il){Rm,9) = 
1pm sin^ 6 and adjust ipm iteratively to give the self-consistent 
dipole moment of the solution, but we encountered numeri- 
cal difficulties with this approach. 

There are two physically plausible choices for the po- 
lar and equatorial boundary conditions: (i) i/)(r, tt) = 
and dip /d9{r,n /2) = 0, as for the initial dipole, or (ii) 
dip / d6{r,±Ti /2) = 0, north-south symmetry. We mostly 
adopt (i) but explore (ii) for completeness in Section 14.61 
where it is shown that field lines on either side of the pole 
are peeled away, leaving the ip = G line isolated, without 
affecting the dipole moment significantly. Strictly speaking, 
the conditions dip /de{r,-K /2) = and dip/dr{Rm,7v/2) = 
force the magnetic field to vanish artificially at {Rn-i,n/2), 
but jm| is affected by less than 0.1% (see Section^J- 



with fi — cos9 and b = — 2r ^(e^ + cot See). Upon solving 
II2II 1 and I22I 1. we obtain 



and 



G{r,n,r',ii') ^^Ng ^ge+i{r,r') 

X (l-M')Cr(M')Cr(M) 



G*{r,tJ.,r',ti') =^Nf ^g*i{r,r') 



(23) 



with 



gi{r,r') = 



(2^-|-l)r'2 ri. 
r< — min(r, r'), r> = max(r, r'). 



g}{r,r') = gi+i{r,r'), 



(24) 



(25) 



(26) 



3 SOLUTION METHODS 

In this section, we discuss three ways to solve 1121 and I14II : 
analytically by Green functions f Section 13.111 . analytically 
in the small- A/a approximation f Section 13.21 . and rmmcri- 
cally, by an iterative algorithm due to iMouschoviaj lll974l) 
(Sections E21 and 0!3J . 



3.1 Green functions 

The Grad-Shafranov boundary value problem 112II . 



AV(n0) = Q(r,e), 



(18) 



with 



ip{R,, 6) =ip, sin 6 and Wm ip{r,e) = Q, (19) 

r — »oo 

can be solved analytically by Green functions if the source 
term Q{r,9) is known as a function of r and 6. In princi- 
ple, Q(r,9) is given by (1121 and 1141 : in practice, it is not 
known analytically. With ip specified on the boundary S of 
the volume V , we can write (see Appendix lAll 



iP{^) = / d^x'G*( 



(20) 



+ / d'yi (ipVG* -G*ViP + hipG*), 
Js 



where G and G* are Green functions for L — fior^ s'm^ dV^ 
and its adjoint L* , satisfying 

d^G , {l-ti^) d^G 1^, 

+ ^2 a. .2 = - ^ - M ) > (21) 



dr 



and 



d^G* ^ (1 - ^^) d^G* ^ AdG' 4^ dG* 

Qj.2 J.2 Q^2 y. Qj. ^2 



(22) 



Ne ^2{£ +!){£ + 2) {2£ + 3)-\ (27) 
and hence, from 1201 . we arrive at the complete solution 

.2\ 



/I roo 
dfi' / dr'r'^geir' ,r)C^^^{fi')Q(r' ,fi'). 
-1 J R, 



(28) 



C^^^lfi) denotes a Gegenbauer polynomial of order £ (see 
Appendix 1X1. 



3.2 Analytic approximation for small Ma 

In the limit of small Ma, where dM/dip and hence Q{r, 6) are 
small, one can show (see Appendix IA31 that the magnetic 
fiux distribution reduces to 



7/'(r,e) = Vi(r,e)(l-&'Ma/Mc) 
far from the star (r 00), with 

Mc = 2TTGM,ipt/fj.ociRl 
For convenience, we write this using CGS units 



Mc 



1.2 X 10" 



108. 



1012G 



(29) 
(30) 

(31) 



where A/* = I.4M0 and R, = let'em. 

It follows that the magnetic dipole moment scales as 
|m| = |mi|(l — Ma/Afc). This scaling agrees, in the small- 
Ma limit, with empirical scalings of the form |m| = |mi|(l -|- 
Ma/Mc ) " \ wi th Mc « lO'^Afp, that have been proposed i n 
the literature dShibazaki et alJll98fll:fGheng fc Zhandll99?^l. 



3.3 Iterative numerical scheme 

To solve and itTHl self-consistently for ip{r, 9) for Ma » 
Afc, we employ an itera tive numerica l algor ithm similar to 
the one introduced by iMouschoviasI (^T^) to study the 
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Parker instability of the Galactic magnetic field. The algo- 
rithm and its performance are discussed in detail in Appen- 
dices^ and |B| and summarized briefly here. 

Given dM/d'4) and an initial guess we cal- 

culate the locations of A'^c contours of ijj, spaced either lin- 
early or logarithmically in capturing topologically discon- 
nected contours and closed loops (|&nvdor 197^. We then 
compute from 1141 and hence F'[^'"'] after poly- 

nomial fitting (simple differencing causes numerical diffi- 
culties; see Appendix IB3^ . The Poisson equation 1)12^ is 
solved with this source term using successive overrelaxation 
to obtain V'ncw('', 6^), and the next iterate is obtained by un- 
derrelaxation: = e^^'-i/)'"' + [1 - e'^'jt/iiol, with 

< 0'"' < 1. Iteration continues until the convergence 
criterion l^/'iow'"'' — ^'"'| < elV'ncw'^'l is satisfied on average 
across the grid. We usually take e = 10~^ in this paper. 
Physically, the algorithm starts from a trial magnetic field 
(and associated current distribution), the accreted mass Ma 
is distributed among the flux tubes according to dM/dij), 
mass is allowed to slide up or down flux tubes to achieve 
hydrostatic equilibrium along B, a new current distribution 
is computed that balances forces perpendicular to B, and 
the process is repeated. 

3.4 Numerical convergence 

There is no general rule for choosing O'"'. We find, 
by experimentation, that one must decrease 1 — 0'"' as 
Ma increases; a useful rule of thumb is 1 — G'"' « 
(Ma/lO-■^M0)-l(V'./lOV'a)"^ for Ma > 10-'^ M0. More de- 
tails can be found in Table IHH and Appendix IB4I At least 
2/[(l — 0)logjo(e)] iterations are required for convergence; 
bootstrapping is recommended, i.e. using the equilibrium so- 
lution for a lower value of Ma as the first iterate instead of 
the dipole. We show in Appendix IB5I that the error in -0 
averaged over the grid, scales as G~^'^ , where G is the num- 
ber of grid cells in each dimension.. The optimum number of 
contours is A'^c ~ G — 1; F'ii^) becomes jagged for Nc ^ G 
due to grid crossings, (as demonstrated in Figure I5T1 . To 
concentrate maximum grid resolution near the stellar surface 
and at the edge of the polar cap (-i/i = i/ia), where screening 
currents predominantly reside and gradients of p and t/j are 
steepest, we scale the r and 9 coordinates logarithmically as 
described in Appendix IB II 

Figure Q displays the mean residual as a function of 
iteration number. Convergence is rapid for Ma < 10"'' Mq 
and poor for Ma > 1O~*M0. Large fluctuations in the mean 
residual are mainly due to the polynomial fit to f (V') and 
the appearance of magnetic bubbles (see Section f4. 711 . 
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Figure 1. Mean residual versus iteration number. Ma = 
IQ-^Mq, e = 0.995 (top). Ma = IQ-^Mq, 6 = 0.95 (bottom). 

spreading of the accreted material and the formation of an 
equatorial magnetic 'tutu', including a criterion for the onset 
of spreading, is discussed in Sectional] The scalings of |ml 
with respect to Afa and b = tpt/i'a. are derived analytically 
and numerically in Sections 14.51 and 14.61 Finally, the forma- 
tion of magnetic bubbles disconnected from the star — a new 
effect — is explored in Section|^7| We start from the undis- 
turbed dipole I16II and adopt the following p hysical parame- 
ters: M, = 1.4Af0, R, = lO^m, B, = 10** T (iHartman et alJ 
ligg?!) . Cs = lOV s'S xo = clRllGMt = 0.54m a nd hence 
a = R^/xo = 1.86 X lO" llBrown fc BildstenllTogjl . The re- 
sults are mostly presented in rectangular (r, 0) plots scaled 
logarithmically where appropriate to emphasize the bound- 
ary layer of compressed magnetic field. 



4 EQUATORWARD HYDROMAGNETIC 
SPREADING 

In this section, we present the results of the analytic and 
numerical calculations described in Section |3 The hydro- 
magnetic structure of the polar 'mountain' formed by the ac- 
creted material is described in Sections 14. Il and l4.3l and com- 
pare d with previous calculati ons in which Flip) is arbitrary 
(e.g. iBrown fc Bildstenlll998ll . The physics of equatorward 



4.1 Structure of the polar mountain and 
equatorial magnetic tutu 

During the early stages of accretion, matter piles up on the 
polar cap, confined by the tension of the polar magnetic flux 
tube. However, for A^a 1O~^M0, the hydrostatic pressure 
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Figure 2. Polar plot of equilibrium magnetic field lines 
(solid curves) and density contours (dashed curves) for Ma = 
10~®Mq. The coordinates measure altitude above the stel- 
lar surface. Density contours are drawn for 7?pmax (maxi- 
mum at the pole, pmax = 2.52 X lO^^kgm"'') with rj = 
0.8, 0.6, 0.4, 0.2, 0.01, 0.001, 10"*, 10"^, lO"'*, lO^^^ 
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Figure 3. Polar plot of magnetic field lines after accretion (solid 
curves) and before accretion (dashed curves) for Ma = IO^^Mq. 
The coordinates measure altitude above the stellar surface. 



at the base of the accretion column overcomes the mag- 
netic tension and matter spreads over the stellar surface 
towards the equator, dragging along frozen-in polar field 
lines. The spreading distorts B, generating screening cur- 
rents fJ-Q^V X B, which act to decrease the magnetic dipole 
moment (|m| is dominated by the polar field). In its turn, 
the spreading is also counteracted by the tension of the mag- 
netic field lines compressed towards the equa tor. These equa- 



tl^ngli 
|2Q0lD, 



greatly 



toria l stresses, neglected in previous work ( Hameurv et alJ 
Il983l: iBrown fc BildstenI ITool iLitwin e 

increase the Ma required to reduce |m|. 

Figure |5| shows the magnetic field configuration and 
density profile for Afa ~ W~ ^Mp, in cross-section (c f. 
schematic version in Figure 1 of lMelatos fc PhinnevlboOll) . 
The 'polar mountain' of accreted material is readily appar- 



ent, traced out by the dashed contours. Figure |3 shows the 
distorted magnetic field configuration overlaid on the field 
lines of the undisturbed dipole. The distorted field exhibits 
a pinched, flaring geometr y, termed an 'equatorial tutu' by 
iMelatos fc Phinnevl J200J) . A more complete view of the 
overall hydromagnetic structure can be gained from Figure 
|1| The tutu-like field is shown again in Figures0Ja) and^Ib) , 
while the polar mountain (p) is shown in Figurelljc). In Fig- 
ure Ud), where the radius of curvature of B is smaller than 
the hydrostatic scale height xq, the toroidal screening cur- 
rents are confined below altitude xq, and are concentrated 
near the polar cap. (Note that ohmic dissipation, neglected 
here, is important at these scales). The J x B force per unit 
volume [Figure 2Ie)] balances the pressure gradient [Figure 
0Jf)], preventing the accreted material from spreading all 
the way to the equator. 

The maximum density, attained at {r,9) — (i?«,0), 
is found empirically to be p, 
10i2(&/10)2(Ma/10-i° 
estimates for p(r, 6) 
ried out in Appendix ^ Consequently, pmax exceeds the 
crustal density 4 X 10^" kgm~^ for Ma > IO'^Mq. In reality 
this overdensity is moderated by sinking iBrown fc BildstenI 
119981: IChoudhuri fc KonaJbOOa) . which is prevented by the 
hard surface in our calculation. (We can alleviate the over- 
density in our model by relaxing the isothermal assumption 
or allowing a nonbarometric density distribution along con- 
tours.) For Afa > IO^^Mq, the maximum magnetic field 
strength becomes unrealistically large (-Bmax 10^^ T) be- 
low an altitude xq, in response to pmax- Suc h field strength s 
formally exceed the yield stress of the crust llRomanMl990l) . 



max - Afa&7(2a^o7ra2) « 6 X 
MQ)kgm~^, in accord with analytic 
Pmaxexp(-a::/2;o) exp(-7/)/V'a) car- 



4.2 Sinking 

The proportion of the accreted material that sinks is not 
well constrained. We consider a crude model of sinking in 
which a proportion s of the accreted matter sinks, leaving a 
proportion (1 — s) to spread. This is modelled by setting 

dM/di, = (Afa/Va) [(1 - s)e-"*/'^V2(l - e-"^- ) + s/{2b)] . 

(32) 

We find that for A'/a = 10"^Mq, for s = 0.9 (i.e. aU 
but 10~^Mq redistributed.), the resultant dipole moment 
increases from « 0.91 to 0.96. 



4.3 Self-consistent mass-flux distribution 

In previous studies of neutron star a ccreti on , F(tl}) was cho- 





1998 
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jw^, ^^^^^^^^ 

Litwin et al. 




20q3 



ar bitraril y dUchida 
il983: .Brown fc B ildstcn 

iMelatos fc Phinnevl (200 ll) . In this paper, by contrast, 
we determine F{tp) self-consistently by solving 1121 and 
1141 1 simultaneously for a physically plausible choice of 
dM/dtp that places most of the accreted material at 
the poles of the initial, undisturbed dipole. Figures (SJa) 
and Ob) compare our self-consistent F{il}) against the 
functional forms guessed by previous authors for two values 
of Ma. The differences are significant, especially near the 
pole (?/) « 0). [An analagous difference was discovered by 
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Figure 4. = 10 ''Mq. (a) Magnetic field lines ('</> contours), (b) magnetic field strength (|B| contours), (c) den- 

sity, (d) current, (c) Lorentz force, and (f) pressure gradients. For each quantity x, values JjXmax are plotted, with rj = 
0.8, 0.6, 0.4, 0.2, 0.01, 0.001, lO"-*, 10"^, lO""^, lO'^^. Maximum values are found to be pmax = l-7xl0"kgm-3, |B|max = 3.9xl0llT, 
|J|max = 2.0 X 10i5Am-2, |J X B|max = 3.3 X lO^^Nm-^, |VP|max = 1.9 X lO^^Nm-s. 
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Figure 5. Magnetic configuration as a function of accreted mass. Ma = 10 ^Mq, pmax = 4.2 X lO-'^'^kgni |B|niax = 3.5 X lO^T 
(top); Ma = 10-'^ Mq, Pmax = 4.2 x IQl^kgrn-^ |B|max = 1.1 X IQl^T (middle); and Ma = 10-^Mq, p^ax = 1.7 X IQi^'kgm-S, 
|B|niax = 3.9 x lO^^T (bottom). Displayed are contours of (/■ (left, solid), |B| (right, solid) and p (left and right, dashed), with values 
rjpmax and ??|B|max, where ?? = 0.8, 0.6, 0.4, 0.2, 0.01, 0.001, lO""*, lO'^, lO''^, lO'^^ 



Burial of the polar magnetic field of an accreting neutron star 9 




tp/Tp Distance from pole, (y [m]) 




Figure 6. Comparison of tlie self-consistent F{ip) (solid) for (a) 
A/a = 2 X 1O~5M0,6 = 10 a nd (b) Ma = 1.5 X IQ-^ M p,,b = 3 
with others in the literature; iBrown fc Bildstenl Jl99d) ( dotted) 
iLitwin et all 120011) (dot-dashed). and lMelatos fc PhinnevI i200l|i 
(dashed). These curves are scaled to the value at the pole, -F(O). 



iMouschoviaj il974) when solving for the final states of the 
Parker instability in t he Galaxy self -consistently, relative 
to the previous guess of iParkeil j^^).] A polynomial fit in 
the numerical code yields the approximate form 



+ 0.21i/i^ -0.021i/i + 0.1333) 



(33) 



for Np = 4,Ma = 1.5 x 10-''A/o,6 = 3, -0 = -i/j/t/^a, and the 
choice of dM/dip in Section f2.2l 



4.4 Onset of spreading 

iBrown fc Bildstenl lll998f) showed that, for an initially verti- 
cal field and neglecting the stress from compressed equa- 
torial flux, the condition for spreading is given by q = 



c; 



0.01, where iJcap = {Ri/R^ 



xl/2 



the polar cap radius. Given that cvmin = -B«/(2/iopmax) ~ 
O.27(B./lO*T)2(Afa/lO"^^M0)-\ it follows that the ac- 
creted matter distorts the magnetic field negligibly for 



Figure 7. Magnetic field lines (solid) and density contours 
(dashed) around the polar cap for Ma = 10~^^Afci and b = 100; 



the parameters used in Figure 2 of lLitwin et alJ 1200 



Afa < 3.7 x lO~''Af0 (i.e. a^in > 0.01). Above this 
value, progressively more distortion takes place for Afa — 
10"^, 10"^, lO'^^Af©, as shown in FigureEl The curved field 
lines have a large tangential component, previously negligi- 
ble, which increases the magnetic field strength — B — sub- 
stantially for Afa > Mc/a = 8 x IO'^^Mq, as predicted by 
the the Green function analysis in Appendix IA3I This in- 
duced magnetic pressure balances the overpressure of the ac- 
creted material. However, the effect on the magnetic dipole 
moment is negligible (< 1%) due to countervailing magnetic 
stresses from the compressed field at the equator; the mag- 
netic radius of curvature is less than xo until Afa exceeds 
1.4 X 10"® Af©, as proved in Appendix IA3I Note that the 
compression of field lines as |B| increases is imperceptible 
in Figure 13 because of the extent of the horizontal axis; |B| 
increases predominantly due to the Bq component. The top 
of the boundary layer is roughly where Bq vanishes, i.e. at 
an altitude x = a;o ln[(a + 1)/(1 - A'/c/Afa)] ~ 5.3 m. for 
A/fa < Mc 

W e zoom in on the pole to compare with 
12001") in Figure Q Our results differ because 
(2001) has a free boundary at the polar cap edge and ig 



Litwin et al 



Litwin et al 



nores the terms in the Grad-Shafranov equation, while 
we impose north-south symmetry at the equator, with no 
condition at the polar cap edge. This allows the compressed 
magnetic field equatorward of the polar cap to push back on 
the polar flux tube. Furthermore, we prescribe dM/dip and 
calculate F{'tp) instea d of prescribing F(^ ). We observe cur- 
vature comparable to lLitwin et al.l 1200 for Afa ~ lO~*Af0 
and thus the ballooning instability may be relevant, but a 
detailed calculation is beyond the scope of this paper. 

An order of magnitude estimate of Afc, the mass re- 
quired to buckle the magnetic field, can be obtained in 
the following way. The hydrostatic pressure at the base 
of the accreted column is given by Pmax = CgPmax ~ 
c^MJP' /{2-KFLixo), i.e. the weight per unit area of the mass 
Afa spread over approximately two hemispheres. This pres- 
sure is balanced by the tension of the magnetic field com- 
pressed into a layer of width L 6 m along the sur- 
face and at the equator (see Figure I^J. In the layer, we 
have Bi ~ B,R,/L ~ 10^^ T by fiux conservation. Hence 



10 D. J. B. Payne & A. Melatos 



i.a [ 



(a) 



the hydrostatic and magnetic pressures balance for Ma = 
2B'iRtxo/{cinoL'^) ~ 2 X 10"'' in accord with the nu- 
merical results. 
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Figure 8. (a) Dipole moment |m| as a function of altitude x = 
r — normalised to |m| at r = R,, for Mh/Mq = 10~^ (top), 
3 X 10-*^, 10-5, 3 X 10-5, 10-4, 3 X 10--* (bottom), (b) Dipole 
moment |m| as a function of accreted mass Ma, normalized to the 
dipole moment |mi| before accretion, for 6 = 3 (crosses) and b = 
10 (triangles), (c) Power law fit (solid), with la errors (dashed), 
to dipole moment as a function of Ma for Afa > 2 X lO^^Afg. 



4.5 Reduction of the magnetic dipole moment 

The magnetic dipole moment 



3»f 

4 



d{cos9) cos 9 Br (r, 



(34) 



is plotted as a function of r in Figure |HIa) . The screening 
currents are confined to a thin layer above the stellar surface; 
|m| is essentially constant with r above this layer. The layer 
is compressed as Ma increases, with half-width comparable 
to xo for Ma « 1O~'*M0. The asymptotic value of |m| also 
decreases as Ma increases, as expected; equatorward hydro- 
magnetic spreading drags magnetic flux away from the pole, 
and |m| is sensitive to Br near the pole through 13411 . 

FigurelSIb) is a plot of |m| versus Ma. For < = 
1.2 X 10"'' Mq, |m| decreases proportional to (1 — Ma^/Mc), 
as predicted analytically in Section 13.21 For Ma > Mc, we 
obtain the empirical relation 



|m|/|mi| = (Afa/4.6 x 10~''Mq)' 



(35) 



by fitting a power law to the numerical results for Afa > 
5 X IO^^Mq, as in Figure [Sl^c'l. However, 13511 cannot be ex- 
trapolated reliably to the regime A'/a > 1O~*'M0 for two rea- 
sons. First, our numerical scheme is limited by the steepness 
of gradients in the source term of I12II . Physically, these scale 
as the hydrostatic pressure, with dF/dip oc M^b as shown in 
Appendix^ We encounter convergence errors above 50% for 
Ma6 > 3 X 1Q-*Mq. Second, it is shown in SectionlTTlthat 
magnetic bubbles, disconnected from the stellar surface, are 
created for Mab > 10~'*AfQ, leading to the steep dependence 
of |m| on Ma in 1351 . A word of caution: when bubbles ap- 
pear, it is unclear how to interpret |m| = |m|star + |m|bubbie. 
In reality, one has |m|bubbic = as r ^ oo, because the flux 
surfaces of the bubble are closed. However, at r = Rm, the 
ingoing and outgoing flux tubes of the bubble do not cancel 
perfectly and one finds |mlbubbic 7^ 0; bubble-related cur- 
rents outside the solution domain (r > Rm), need to be 
included in order to recover I ml bubble ~ 0. 



4.6 Polar cap radius 

We now discuss the effect on |m| of varying b — tpa/tp*, or 
equivalently the polar cap radius Reap = R, sm~^(b~^^^). 
Although 7?cap is not known exactly without a detailed 
model of the fl ow of matter from the accretion disc to the 
stellar surface iArons et al ]ll984'), estimates of its size from 
1151 1 are typically of orde r 1km [.Litwin ct al. 2001 j, i.e. 
b > 100, for = 10* T llArons et alj(l984ll . Figure H il- 
lustrates the magnetic configuration obtained for 6 = 3 and 
b — 10, and |m|/|mi| is plotted versus Ma in Figure |HIb), 
denoting fe = 3 by crosses and 6 = 10 by triangles. Note how 
the equilibrium state changes with b. The mass-flux distri- 
bution dM/dtp oc exp(— i/;/i/)a) implies a surface pressure dis- 
tribution F{tp) oc bexp{—ip/tpa), so larger b means steeper 
pressure gradients. (Numerical difficulties set in for b = 30, 
which can be partly alleviated by stretching the coordinates 
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around ipa.) Importantly, however, we find Mc and the dipole 
moment |m| are independent of b. This is supported by the 
Green fijnction analysis in Appendix^ except when dM/dif> 
depends explicitly on b, e.g. dM/dip oc 6(1 — ip/bipi^)^^^ im- 
plies |m| oc 6~^. 

A related issue is whether |m| is affected in an unre- 
alistic way by the boundary condition ?/) = at S = in 
the examples presented so far. It is conceivable, for exam- 
ple, that the -0 = line is unstable ('on a knife edge'), while 
neighbouring field lines are peeled away by accretion, unless 
it is forced to remain rectilinear artificially. As it happens, 
however, this is not the case. Figure [Till shows the output of 
an experiment where the grid extends from equator to equa- 
tor (1^1 < 7r/2), and no boundary condition is imposed at 
6 = 0. Clearly, the magnetic field and density profiles remain 
symmetric about the magnetic pole, with ?/j = at 6 = 
emerging naturally, while |m| is essentially unchanged. 



4.7 Buoyant magnetic bubbles 

From the sequence of panels in Figure |S] (Ma/M© = 
10"^ 10-^ 10"^), we observe that the magnetic field be- 
comes increasingly distorted as Afa increases. Eventually, for 



Afa ^ W~^Mq, closed magnetic bubbles are created that are 
disconnected topologically from the surface of the star. This 
phenomenon is illustrated in Figure fTTI At the value of Ma 
where a bubble is first created, a magnetic neutral point (Y 
point) is observed to form on a field line near, but not at 
the pole. The bubble closes at r < Rm in our simulation, 
but this may be a result of the approximate free boundary 
condition dip /d9{Rm,d) = 0; in reality, it may connect to 
the accretion disc. 

Bubbles correspond to a loss of equilibrium, analo- 
gous to that which occurs du ring eruptive solar phenomena 
jKlimchuk fc Sturrocld Il989^ ■ where no simply connected 
hydromagnetic equilibrium exists. In the Grad-Shafranov 
boundary problem, the source term oc F'{ij}) in in- 
creases with MJ}, boosting A'^i/) and hence ip above the sur- 
face (to balance the weight of the added material through 
the Lorentz force). Above a critical value of Ma, fiux surfaces 
are created with ip < Q ox ip > ip,, which are disconnected 
from the star and either form closed loops or are anchored 
'at infinity' (here the accretion disc). This is shown explic- 
itly by IIA19II in the special case F'{ip) = constant. From 
our numerical results, we conclude that the critical Ma for 
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Figure 10. Equilibrium magnetic configuration in the north- 
ern hemisphere, showing contours of constant ip {solid) and p 
(dashed). Equations 1121 and 1141 are solved in the domain 
\9\ < tt/2 here, compared to < S < it/2 in earlier figures, in 
order to test the validity of the ^ = boundary condition at 
9 = 0. 

bubble creation satisfies 

Ma > 10"''&"'Mq . (36) 

Note that bubble creation is a topological imperative. It 
is not the result of a hydro maRnetic instability e.g. in- 
terchange or Rayleigh- T aylor iBernstein et aLlFlOBSl : iParkerl 
Il966l: lMouschoviaslll974l) . 

Are the bubbles merely numerical artefacts 
jBrown fc BildstenI Il99d) ? No. Equation JKwi demon- 
strates explicitly that flux surfaces with i/; < or i/i > 
are created for Ma satisfying 13611 . at least for F'{tp) = 
constant. A more subtle issue is whether bubbles are the 
by-product of an artificial assumption in our idealized 
calculation. For example, if submergence of accreted 
material were permitted, it might reduce the pressure 
gradients that produce the bubbles; on the other hand, 
ohmic dissipation would facilitate detachment of bubbles in 
a pinched, Y-point configuration. 

On some runs, bubbles appear and disappear during 
the iteration process. This happens because the mass-flux 
distribution is not conserved inside a bubble, although the 
code attempts to maintain flux freezing at the edge. If the 
route to convergence is a rough proxy for ti me-dependent 
behaviour, as argued by iMouschoviasI il974ft for iterative 
relaxation algorithms, the appearance and disappearance of 
bubbles may represent evidence — though not proof — of 
transient evolution in reality. 

As the bubbles are disconnected topologically from flux 
surfaces anchored to the star and accretion disc, they do 
not contain any accreted material (in the ideal-MHD limit 
of zero cross-field transport) and are lighter than their sur- 
roundings. It is therefore possible that they rise buoyantly 
and ultimately escape the magnetosphere of the neutron 
star. This possibility cannot be investigated rigorously in 
the context of the equilibrium calculations in this paper; it 
is considered qualitatively in Section [5.31 
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Figure 11. Magnetic configuration for Ma = 2 X 10~^Mq, 
b = 10, showing the creation of a bubble. We plot contours 
of constant (top) for initial (dashed) and final (solid) states, 
and final |B| contours r7|B|inax (bottom, solid) and p contours 
(dashed) 7?Pmax; Pmax = 5 X 10l''kgm-3, Bmax = 6.3 X 10" T, 
and -q = 0.8, 0.6, 0.4, 0.2, 0.01, 0.001, lO"*, 10"^, 10"'^, lO-^^. 



5 TIME-DEPENDENT EFFECTS 

In this section, we discuss critically (but qualitatively) how 
the results of this paper may be affected by time-dependent 
processes that cannot be modelled by a quasi-static sequence 
of hydromagnetic equilibria. We consider Parker instabili- 
ties in Section 5.1, ohmic dissipation in Section 5.2, and the 
buoyant rise of magnetic bubbles in Section 5.3. 

5.1 Hydromagnetic instabilities 

The computed equilibria are manifestly distorted. Buoy- 
ancy of the compressed magnetic fiux can drive long-wave, 
slow MHD modes that overturn the accreted matter on 
the Alfven time-scale ta, as in the global R ayleigh- Taylor 
instability of the Galactic magnetic field jParkeJ 1 19661 : 
IMouschoviasI Il973) . When the accreted matter bends the 
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polar magnetic field, there exists a significant component of 
magnetic field perpendicular to gravity, a condition for the 
onset of the Parker instability. In a plane-parallel geometry, 
wavelengths longer than A — 4:TTxo/{2a + 1) are unstable 
iMouschoviaslllQT^l . where a = /{2^op). The geometry 
of an accreting neutron star is far from plane-parallel. Never- 
theless, in hydromagnetic equilibrium, one has a ~ 1 locally 
in the boundary layer, yielding A = 2 m. The failure to 
converge at large Ma is also a hint, though not a proof, that 
the Parker instability may operate. iMouschoviaj il974l) ad- 
vanced a similar convergence-based argument for the stabil- 
ity of a stratified vertical column with periodic boundaries. 

5.2 Ohmic dissipation 

Our calculations are performed under the assumption of 
infinite conductivity and hence flux freezing. In reality, 
the accreted matter is resistive due to electron-phonon 
and electron- impurity scattering jBrown fc BildstenI ITool 
ICumming et al. 2001), potenti ally enhanced by accretion- 
induced heating (Romani 1990: ^ Urpin fc Geppertiri995ll . 

The ohmic dissipation time-scale for a fiux tube of width 
L is given by Td = /locrL^, where a is the electrical conductiv- 
ity. For typical conditions, we take a ~ lO^^s"^ and hence 
obtain Td ~ 10^^(I//ii, )^ s. In comparison, the flow time- 
scale is given by Tf = AuR^pL / M^., where Ma is the accretion 
rate. For p = 4 x 10^*kgm~^, L = 6 m (Section and 
Ma = 1 X lO"*Af0 yr"\ we find rt = 1 x 10^ yr. Therefore, 
for the length-scales characteristic of the compressed flux 
layer, we have t{ > Td and magnetic flux diffuses through 
the accreted material, broadening the compressed flux layer 
until it is thick enough (L ~ 600 m) that Td t/ and 
further thickening ceases. iBrown fc Bildsteiil (|W9^ showed 
that Tf/Td depends only on Ma and not on depth in the 
crust. 

Note that the buried flux is resurrected on the time- 
scale Td, after accretion stops. As the noted dipole field re- 
asserts itself, |m| increases. A full analysis of this process is 
left to future work. 

5.3 Buoyant bubbles 

Closed magnetic bubbles have p = within and are 
lighter than their surroundings, as discussed in Section f4. 71 
They tend to rise buoyantly at the local Alfven speed 
VA = (-B^/pop)^'^^ = Ca(2a)^'^^, and hence escape the 
magnetosphere in « 2 yr. Assuming an accretion rate 
of lO"*M0yr-\ it takes ^ lO''6~^(Afa/lO"*M0 yr"^)"^ 
yr to accrete enough mass to create one bubble. More- 
over, a typical bubble encloses ~ -i/ja of magnetic flux. 
Hence we conclude that magnetic flux is being expelled 
episodically from the magnetosphere at an average rate of 
10"(V'*/10^'^Tm2) (Ma/lO-*M0yr-i)-i T m^ yr'V Note 
that expelled flux is subsequently replenished by the current 
deep in the star (since i/) is flxed at the surface). 



6 CONCLUSION 

Observations of low-field binary neutron stars and recycled 
pulsars imply that the magnetic dipole moment of a neutron 
star is reduced by accretion. In this paper, we undertake 



a self-consistent analysis (numerical and analytic) of one 
mechanism that may account for the reduction observed: 
polar magnetic burial in the ideal-magnetohydrodynamic 
regime. Our analysis has several new features, (i) Flux freez- 
ing is strictly enforced when connecting the final and ini- 
tial magnetic configurations by solving self-consistently for 
the mass-flux distribution rather than specifying it ad hoc 
(see Figure ISJ. (ii) The Lorentz force due to equatorial mag- 
netic field lines compressed by equatorward hydromagnetic 
spreading is included when calculating the confinement of 
the polar accreted column, (iii) Numerical methods are de- 
veloped for treating accreted masses up to « 1O~*M0, (cf. 
10"^" A/0 in previous work), where the field is dramatically 
distorted and high-order multipoles dominate. 

We report two key results, (i) Afa ^ 10~^Mq must 
be accreted in order to reduce significantly the magnetic 
dipole moment |m| of the star, contrary to previous esti- 
mates (Ma ~ 10"^" M©) which neglected equatorial mag- 
netic stresses. For small Ma, we find |m| = |mi|(l — Ma /Mc), 
with Afc = 1.2 X 1O"^M0, (cf. IShibazaki et al.lll989l) . For 
A^c < Ma < 1O-*M0, we find |ni| = |mi| (Ma/4.6 x 
1O~^M0)^'^^*'''^^ . (ii) When enough mass is accreted, such 
that Ma > 1O~*&~^M0, the hydrostatic pressure gradient 
generates flux surfaces with ?/) < or i/> > i/)*, creating 
closed magnetic bubbles that are disconnected topologically 
from the stellar surface. The bubbles are valid solutions of 
the Grad-Shafranov boundary problem, as confirmed by an- 
alytic. Green function calculations; they are not numerical 
artefacts or fingerprints of hydromagnetic instabilities (e.g. 
Parker) . 

Several of our assumptions need to be relaxed in future 
work, including (i) perfect conductivity, (ii) an impenetra- 
ble stellar surface, and (iii) axisymmetry (which tends to 
suppress hydromagnetic instabilities, for example). Finally, 
the uniqueness of the hydromagnetic equilibria we compute 
numerically is yet to be established. 

Our equilibria serve as useful starting points for explor- 
ing the stability of the magnetic configuration during and 
after accretion. Our theoretical results will be tested against 
observational data from binary neutron stars and recycled 
pulsars in a companion paper. 
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APPENDIX A: ANALYTIC SOLUTION OF THE 
GRAD-SHAFRANOV PROBLEM 

In this appendix, we solve the Grad-Shafranov equation 
II12II . together with the boundary conditions II17II . by a Green 
function approach. 



Al Green theorem for the Grad-Shafranov 
operator 

An operator Ltli = V^ip + b ■ Vip + ap, acting on a function 
ip, possesses an adjoint L*ip = \/^ip — V • (hrp) + cip- The 
Grad-Shafranov operator L — figr^ siri^ OA^ , defined by Q 
in spherical polar coordinates, has b — —2r~^{er + cot 9eg) 
and is not self-adjoint (cf. V^). Letting G and G* be the 
Green functions associated with the operators L and L* re- 
spectively. 



LG'(x,x') = (5(x-x') 



L*G* {x., x') = (5(x — x' 



(Al) 



(A2) 



related by the reciprocity relation G*(x, x') = G'(x',x), we 
arrive at the Lagrange identity 

iPL'-G* ~ G'LiP = V ■ {ipVG* - G'ViP + hiPG"). (A3) 

Upon integrating l|A3fl over a volume V , bounded by a sur- 
face S, and using the divergence theorem, we obtain 

(V'L*G* - G*L'4))dV 



{ipyG' - G'VV- + hipG*) ■ ndS, 



(A4) 



where n is the unit vector normal to S. Given a boundary 
value problem Ltp{x) — Q(x) in the volume V, with tp given 
on the boundary 5", we combine <A2t and IIA4ll to obtain 



V(x) 



d^x'G'C 



+ [ d^x'n-(V'VG*--G*VV' + bV'G*). 
Js 



(A5) 



A2 Green function for the Grad-Shafranov 
equation 

We wish to solve 1121 for tp in r > subject to Dirichlet 
boundary conditions (1171 on tp a,t r — and r — > oo. In 
cylindrical symmetry, the volume and surface integrals in 
1201 1 reduce to surface and line integrals respectively. It is 
convenient to make the substitution /i — cos9, whereupon 
H12|l becomes 



+ 



nor (1 - A* ) 



dF{iP) 
dip 



(A6) 



X exp(-(^o/Cs - GMr/Rtc^) 



and dF/dxp is a function of r and through ip{r,jj,). 

We redefine L to be the operator on the left-hand side 
of 1IA6I 1. and Q{r,iJ,) to be the source term on the right- 
hand side, known explicitly once F{'tp) is known. The Green 
function G for L satisfies 



d^G il-fi'')d^G 



+ 



^ 5{r-r')S{n-n') (A7) 



and the Green function G* for L* satisfies 

a^G* {l-fi^)d^G'' AdG* 4fidG* 



g^2 



r dr i9/i 



(A8) 



Equation IA8II is separable. We expand the solution in terms 
of orthogonal Gegenbauer polynomials G|^^(/i), viz. 

00 

G{r,^l,r' ,11') = y^gi{r,r') 

7^0 (A9) 



with 



d^ge{r,r') £{£+!) , _2., rATn\ 

:h — -" 2 9e[r,r)=r S(r-r). (AlO) 

dr^ r'^ 



The Gegenbauer polynomials satisfy 



+ £(£ + i)[(i-m')g1/1(m)] = o 



(All) 
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and are related to associated Legendre polynomials via 



20/2^3/2 



C;^ii(^i). The first few are listed for ref- 



erence: Co^\^l) = 1, Cf^ifi) = 3m, ^(^i) = 1(5^^ - 1), 



Cl^'M = I (7m' - 3m), Cr^(M) = f (21m' - IV + 1), 



-,3/2 



C5 (m) = |(231m - 210^'^ +35At). They satisfy an orthog- 
onality condition: 



with 



y (1 - M')C^'"(M)C^V'(/i)dM = A^^-Jf. 

2(£ + l)(^ + 2) 



(21 + 3) 



(A12) 



(A13) 



We solve lAlOt for ge{r,r') subject to the follow- 
ing conditions: (i) gi{r,r') is continuous at r = r', (ii) 

= r'-2, (iii) ge{R„r') = 0, and (iv) 



hm [ 

^0 L 



d3l(r,r') 



lim gi{r,r') — 0. The result is 



e{r, r') 



21+1 



(2£-H)r'2 ri 
with r< — min(r, r') and r> — max(r, r'), yielding 

X (l-M')Cr(M')Cr(p)- 

Similarly, we obtain 

oo 

G*{r,ii,r' =^N;r^g*e{r,r') 

x(l-M")cf^M')cf'(/^), 



(A14) 



(A15) 



(A16) 



with ge{r,r') — { — )'^gi+i{r,r'). Equations HA 1511 and 
JA16I I are consistent with the reciprocity relation. Note 
that basis functions for G and G* for a non-self-adjoint 
operator are mutually but not individually orthogonal 
jMorse fc Feshbachlll953l) . 

Upon combining the Green theorem ljA4^ , the definition 
of G* <A16ll . the boundary conditions i/)(i?*,M) ~ tp*{l — 
m'), tp{r,±l) = 0, lim t/.(r,M) = 0, G* (R,, fi,r' , fi') = 0, 

r — >oo 

and the surface gradient Vil^{r, ±1) -e^ = 0, we find that the 
boundary integral over C reduces to J^^tp{Rt, n)^^dfj,' , 
yielding the complete solution 



(1-M^) 



X / dfi' 

1-1 J R, 



(A17) 

dr'r'^gi{r' ,r)C^^^ {lJ,')Q{r' , n'). 



A3 Small-A/a limit: Constant Source Term 



To explore the form of the general solution H28|l , in the small- 
Ma limit, we linearise Fi^l)). We consider the special case 
F{^) = Qoi^p, - ^), giving 0(r,M) = Qo{l - M')r'e-'-. By 
orthogonality, only the £ — term survives. We find 



/ '\ /4 —r' J I 

gi(r, r )r e ar 



-[/i(r)-/i(i?,)], (A18) 



with /i(r-) = (r^ -I- 4r^ + 8r + 8)e and hence 



(A19) 

using A^i = 4/3. It is immediately clear that negative values 
of ip are possible if Qo > 4'i/'*-R* [/i (r) — /i(_R,)] raising the 
possibility of closed magnetic loops (bubbles) constructed 
from fiux surfaces tp < or > tl), and hence not anchored 
to the stellar surface. 

In dimensionless coordinates, setting F{^p) = k{b — ip) 
we arrive at 



^(i,M) = ■</'i(i,M) {1 + 



b (A20) 

X [/i(i)e-*-/2(a)(l~e-*)]}, 

where fi{x) = 3m"^ -I- a"^(3i^ -I- 85) -I- a"^(i^ -f 4i^ + 8£) 
and /2(a) = 1-1- 4a"^ -I- 8a"^ -I- 8a"^. For £ < a, where 
the screening currents dominate, 1/) reduces to tpi{x,p,) [1 — 
fcQoa^&"^(l - e"*)]. In Appendix 1X1 we show that k « 
h/{2-Ka?), so the reduction factor is (1 — b^M^/Mc), remem- 
bering that Qo oc Ma. For neutron star parameters, one has 
Mc ~ 1.2 X 1O"*M0. 

We can estimate the thickness of the compressed fiux 
layer, Xh, in the small- Ma regime by solving ip{xh,fl) ~ 
to give Xi, = — ln(l — Mc/Ma) We may also estimate when 
B| changes significantly from its initial value |Bi|. Near the 
surface, the principal component 



{Bi)^.[ab'A'h/M,e-^ + 1 - b^M^/M,{l - e"")] (A21) 



increases significantly for M^/a > 8 x 1O~^^M0, consis- 
tent the numerical results in Section |21 Setting = 
gives an alternative estimate of the altitude below which 
the screening currents are confined, with x = ln[(a -)- 1)/(1 — 
&^ Mc/Ma)] =9.8 consistent with lIMB . Including also the 
radial component Br, we obtain 



B=Bi[4M'(l-&'Ma/Mc(l 



If 



2,1/2 



+ (1 - M )(1 + ab'Nh/M,e-y] 
which reduces to 

B = Bi[An^ + (1 - i^'){ab''M^/M,f]^''^ 



(A22) 



(A23) 

near the surface for 10"^" < Ma/Mc < 10^*^. There is also 
a boundary layer at the magnetic pole, where increases 
rapidly from zero over a short distance. The width of this 
polar boundary layer may be estimated by setting B^ ~ Br, 
yielding r; ttTI. tan-^[2Mc/(a6^Ma)] ~ 20 m. 



A4 Dipole field 

A useful analytic approximation to the source term dF/dip 
can be derived for the dipole field in the early stages of 
accretion. In dimensionless coordinates, defined in appendix 
ini lltil becomes 



il)i{x,jl) ^b{l - + x/a) ^ 



(A24) 



with a = R*/xo and b — tp,/tpa_. From 1141 . we write 
dM/di! = F{i))I{'4)), with 



/(V-) = 27r 



/ dS(l -/i^)'/'(i + a)e-*|V7/ir' 
Jc 



(A25) 
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and 



OjO 
O 



0.5 



0.0 



-0.5 



-1.0 



-1.5 




-5 -4 -3 -2 -1 
log.o (1 - f/b) 

Figure Al. The function Jitp) (solid) compared with J{ip) cal- 
culated numerically for Ma = 10~^AfQ, 6 = 10 (dashed). 



C is a contour of constant ^p, along which we may write 
/i = [1 — '^r/{ab)]^^^; the integral terminates at i = on 
the surface and x — a{b/tp — 1) above the equator. Upon 
rearranging we obtain 



I{i>) =-naJ{i,)b'^^'^{b-'4>) 



-1/2 



(A26) 



with 



Jill;) = / dx{l + x/af 

Jo 



a{b - ip) 



-1/2 



(A27) 

The function J{tp) is plotted in Figure Ell We observe that 
J ~ 1 for all ^ except near the equator. Its limiting be- 
haviour is /(i/i) ■Ko?b~^ as ■4) ^Q, and I{Tp) 2™^ (6 - 
^-)i/2^-3/2 ^ ^ If we choose dM/d^ = exp{-i!)/2 
specifically, then from <A25t we have 

= TTT? exp(-^)(l - ^/by^^lJW]-'. 



27ra2 

Upon differentiating with respect to tp, we obtain 

dF be-''' 
d^ 



(A28) 



27ra2 

X [l + (l-2^)/(2b)][J(V;)]-' 
- {l-i,/b)'^^[JWr''j' 



(A29) 



APPENDIX B: ITERATIVE NUMERICAL 
SCHEME 

Bl Dimensionless equations and logarithmic 
coordinates 

It is convenient to convert to dimensionless variables x = 
(r - R*)/xo, = VM, M = M/Ma, jl = cos6,F = F/Fo 
and B — B/Bo where a;o = c^Ri/{GM,) is the pressure 
scale height, a — R,/xq, Fq — Mi^c^/xq, and Bo — V'a/''o = 
a^B,/2b. Equations 11211 and I14I I take the forms 



^(l-/i')(i + a)2e-*^ (Bl) 
dtp 



dM 



^27TFii,) / ds(l-/i^) 
dtp Jc 



2^1/2 



{i + a)e ^^IVipl' 



(B2) 



respectively, with Qo = imqXqM^c^ /tpl^. Note that /iq denotes 
the permeability of free space (SI units). 

As A/a increases, a thin boundary layer of screening 
currents form s near the surface of the neutron star (see 
Sectional (Melatos fc PhinnevllioOlh . To concentrate max- 
imum grid resolution at the boundary layer and at the edge 
of the polar cap {%p = tps), where the gradients of p and tp are 
steepest, we scale the r and 6 coordinates logarithmically, 
such that 



xi = log(5 4- e + Lj: 



yi = -iog[i 



(l-e- 



(B3) 



(B4) 



must be chosen sufficiently small to ensure at least sev- 
eral nodes per hydrostatic (or hydromagnetic) scale height. 
To resolve steep gradients at the equator, a similar transfor- 
mation exists. 



B2 Grid & Poisson Equation 

We use a grid of {Gx,Gy) cells in {r,9) and Nc contours 
of tp, and choose Nc < Gx to avoid zig-zags in I{tp) due 
to grid crossings which damage convergence, (Figure inU- 
We typically set Gx = Gy — 256 and Nc — 255, with the 
contour values chosen to lie between grid points at the stellar 
surface. 

Given Ma and hence dM/dtp, and starting with a guess 
tp'-°\r,e), the left-hand side of (EU. A^-f/)'"', can be cal- 
culated at each grid point. The resulting Poisson equation 
is solved using a success ive over-relaxation procedure with 
Chebychev acceleration jPress et alJll992|) . stopping when 
the mean residual over the grid satisfies (Atp/tp) < e = 10~^. 

B3 Grad-Shafranov source term and contouring 

To find the source term on the right-hand side of 112II . F{tp) 
is calculated from 1141 . The integral along tp contours re- 
lies on a contouring algorithm adapted from iSnvd er (197^), 
which can follow closed loops and topologically disconnected 
contours. Numerical difi'erentiation of F{tp) by first or sec- 
ond order differencing leads to numerical problems, magni- 
fying small fluctuations in I{tp) and hence F{tp). We over- 
come this by smoothing F{tp) at each iteration step, fitting 
an order Np polynomial (TVp — 10 typically) to [/(i/))]~^, 
viz., [litp)]"^ = X]i=o '^'(I'^S V")* S'lid then differentiating 
F{tP) = [I{tP)]-'-dM/dtP analytically. 

B4 Underrelaxation 

The contour values of dF/dtp are mapped onto the grid by 
linear interpolation and fed into 112t by under-relaxation. 



exph(0-^)/c,^], 



tP 



(n + l) 



e(")^(") + [1 _ e^">]tpi:^ 



ICW ) 



(B5) 



(B6) 
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b 




M^/Mq 






10-10 
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5 10-* 


3 


0.0 


0.5 0.9 


0.995 


10 


0.0 


0.8 0.99 


0.99999 


30 


0.0 


0.99 0.999 





Table Bl. Optimal e'"' values as a lunction ol Ma and b. 



Grid size 


mcan(Ai/)/t/>) 


max{Ai/)/i/>) 


8 


0.0603082 


0.229865 


16 


0.0166787 


0.121939 


32 


0.00307997 


0.0222018 


64 


0.00129922 


0.0328857 


128 


0.000404714 


0.00589775 


256 


0.000286883 


0.0104630 



Table B2. Average and maximum errors as a function of grid 
size. 



is satisfied for all grid-points. As a rule of thumb, we take e = 
; our solutions do not converge reliably 
for Ma > IQ-^Mq. The stability of the solution is checked 
by perturbing it slightly and looking for re-convergence, or 
by resetting 9*"' « 0. 

The optimal value of 0'"' is governed by maximum ip 
gradient between grid-points, which in turn depends on Ma. 
We adjust 9'"' towards unity when IV^new^' - tends 
to increase. Table [54l shows approximate optimal values of 
for different input parameters, chosen to minimize the 
number of iterations while still achieving convergence. 

B5 Testing and convergence 

We tested the code by successfully reproducing the final 
equilibr ium states of the P arker instability, plotted in Fig- 
ure 2 of lMouschoviaj il974) . to an accuracy of 1%. We also 
tested the code against the exact analytic solution lAlQt in 
Appendix^for constant dF/dtp. Table [521 shows a compar- 
ison of the mean and maximum errors as a function of grid 
size, relative to 1A19II . when solving the Poisson equation 
directly [9*"' = 0]. 



Figure Bl. Comparison of convergence for = 0.995, 6 = 3, 

A/a = 8.5 X 10-^Mq- (i) G = 64, Nc = 63 (top), (ii) G = 64, 
Nc = 255 (middle), (iii) G = 256, Nc = 255 (bottom). 



where ■i/'ncw is a provisional iterate and < 9^"' < 1 is 
the under-relaxation parameter at the ri'^ iteration. Con- 
vergence is reached when 

i^i:+'^-^'"^i<ei^i:+''i (B7) 



